$NOFLOATCALLS
$NODEBUG
$STORAGE:2
c****  User's calling program
c****  NEQN=4 so dimension of work array = 4*17 = 68

       implicit double precision (a-h,o-z)
       dimension y(4),work(68),icom(4)
       external orbit
       common alfasq

       open(2,file=' ',status='new')
       icom(1)=0
       write(*,*) 'imeth=,tola=,tolr='
       read(*,*) imeth,tola,tolr
       ecc=0.25d0
       alfa=3.141592653589d0/4.d0
       alfasq=alfa*alfa
       neqn=4
       hstart=0.01d0
       y(1)=1.d0-ecc
       y(2)=0.d0
       y(3)=0.d0
       y(4)=alfa*dsqrt((1.d0+ecc)/(1.d0-ecc))
       x0=0.d0
       xb=0.d0
       icom(2)=0
       icom(3)=0
       do 20 j=1,24
	  xa=xb
	  xb=0.5d0*dble(j)+x0
	  write(2,100)xa,y(1),y(2)
	  call runkut(xa,y,xb,neqn,tola,tolr,hstart,work,
     &		       imeth,ierror,icom,orbit)
	  if(ierror.GT.1) then
	    write(2,100)xb,y(1),y(2)
	    write(2,*) ' ERROR-Problem too stiff or is discontinous'
	    close(2)
	    stop
	  end if
20     continue
       if(icom(4).GT.0)write(2,*) ' Severe round-off error possible'
       stop
100    format(F5.1,3F15.9)
       end
c********************************************************************
c****  User supplied subroutine that contains the system of
c****  differential equations to be solved.
c****  Notice that in this routine it is necessary to have an array
c****  yprime(neqn)

       subroutine orbit (x,y,yprime,neqn)
       implicit double precision (a-h,o-z)
       dimension y(neqn),yprime(neqn)
       common alfasq

       r=y(1)*y(1)+y(2)*y(2)
       r=r*dsqrt(r)/alfasq
       yprime(1)=y(3)
       yprime(2)=y(4)
       yprime(3)=-y(1)/r
       yprime(4)=-y(2)/r
       return
       end
